home *** CD-ROM | disk | FTP | other *** search
/ Practical Algorithms for Image Analysis / Practical Algorithms for Image Analysis.iso / TARFILE.GZ / tarfile / ch_5.4 / fitcrit / wwtheta.c < prev    next >
Encoding:
C/C++ Source or Header  |  1999-09-11  |  11.0 KB  |  342 lines
  1. /* WWTHETA:     function calculates theta plot, finds ww features on it
  2.  *            and returns number of features found
  3.  *                  usage: nThetaFeats = wwtheta (data, edge, wwPar, 
  4.  *                                      nodesInLine, nNodesInLine, nodeList,
  5.  *                                              cycleFlag)
  6.  *              Note that a flag value of BORDERTHETA is put at the positions
  7.  *              at either end of the theta array to indicate the border
  8.  *              theta values which are not valid values.
  9.  */
  10.  
  11. #include <stdio.h>
  12. #include <stdlib.h>
  13. #include "ww.h"                 /* wing-walk include file */
  14.  
  15. extern long thetafeat (struct point *, struct edge, struct theta *,
  16.                        long, struct wwPar, short);
  17. int thetaarc (struct point *, struct edge, struct theta *, long);
  18. double thetacalc (struct point *, long, long, long, long);
  19. long thetasmth (struct theta *, long, long);
  20. double median (struct theta *, long);
  21. long wwends (struct point *, struct edge, struct theta *,
  22.              long, long, long *, long *);
  23. long
  24. wwtheta (data, edge, wwPar, curveFlag)
  25.      struct point *data;        /* sequence of coordinates of lines */
  26.      struct edge edge;          /* lower/upper indices of segment in data */
  27.      struct wwPar wwPar;        /* wing-walk parameters */
  28.      short curveFlag;           /* =0 for curve fit by 2 lines; =1 for curve */
  29. {
  30.   struct theta *theta;          /* theta plot curv. and arc length array */
  31.   double thetacalc ();
  32.   long dataMid,                 /* middle pt on curve */
  33.     iLower, iLow,               /* data indices of lower W-segment */
  34.     iHigh, iHigher,             /* data indices of higher W-segment */
  35.     nTheta,                     /* no. values in theta plot */
  36.     iTheta,                     /* increment in theta array */
  37.     flag,                       /* tests for beginning/end of data curve */
  38.     nThetaFeats,                /* no. of features found on theta plot */
  39.     i;
  40.  
  41. /* if line is too short to fit ww to, forget it */
  42.   nTheta = edge.iHigh - edge.iLow + 1;
  43.   if (nTheta <= wwPar.l)
  44.     return (0);
  45.  
  46. /* allocate space for theta plot */
  47.   if ((theta = (struct theta *)
  48.        calloc (nTheta, sizeof (struct theta))) == NULL) {
  49.     printf ("CALLOC: not enough memory -- sorry");
  50.     return (-1);
  51.   }
  52.  
  53. /* set theta to border flag, and calculate arc length along theta plot */
  54.   for (i = 0; i < nTheta; i++)
  55.     theta[i].theta = BORDERTHETA;
  56.   thetaarc (data, edge, theta, nTheta);
  57.  
  58. /* find start point where WMW fits within beginning of curve */
  59.   for (dataMid = edge.iLow + 1, iTheta = 1;; dataMid++)
  60.     if (wwends (data, edge, theta, dataMid, wwPar.l, &iLower, &iHigher)
  61.         != -1)
  62.       break;
  63.  
  64. /* calculate theta plot */
  65.   flag = 1;
  66.   iTheta = dataMid - edge.iLow;
  67.   while (flag >= 0) {
  68.     wwends (data, edge, theta, dataMid, wwPar.m, &iLow, &iHigh);
  69.     theta[iTheta].theta = thetacalc (data, iLower, iLow, iHigh, iHigher);
  70.     dataMid++;
  71.     iTheta++;
  72.     flag = wwends (data, edge, theta, dataMid, wwPar.l,
  73.                    &iLower, &iHigher);
  74.   }
  75.  
  76.   thetasmth (theta, nTheta, wwPar.m);
  77.  
  78. /* determine WW features from theta plot */
  79.  
  80.   nThetaFeats = thetafeat (data, edge, theta, nTheta, wwPar, curveFlag);
  81.  
  82.   free (theta);
  83.   return (nThetaFeats);
  84. }
  85.  
  86.  
  87. /* THETAARC:    function computers arc length for each theta plot value
  88.  *            computed from the first data coordinate
  89.  *                      usage: thetaarc (data, edge, theta, nTheta)
  90.  */
  91.  
  92. int
  93. thetaarc (data, edge, theta, nTheta)
  94.      struct point *data;        /* data coordinate array */
  95.      struct edge edge;          /* lower and upper <data> indices */
  96.      struct theta *theta;       /* theta plot with arc lengths */
  97.      long nTheta;               /* no. values in theta plot */
  98. {
  99.   double incr;                  /* increment = 1 (horiz/vert), = 2 (diag) */
  100.   long i, j;
  101.   struct point diff;            /* x and y differnces */
  102.  
  103.   theta[0].length = 0.0;
  104.   for (i = 1, j = edge.iLow + 1; i < nTheta; i++, j++) {
  105.     diff.x = data[j].x - data[j - 1].x;
  106.     diff.y = data[j].y - data[j - 1].y;
  107.     incr = (double) (ABS (diff.x) + ABS (diff.y));
  108.     incr = (incr == 1.0) ? 1.0 : DIAG;
  109.     theta[i].length = theta[i - 1].length + incr;
  110.   }
  111.   return (0);
  112. }
  113.  
  114. /* THETACALC:   function calculates theta angle between two line segments
  115.  *                    usage: theta = thetacalc (data, iLower, iLow, 
  116.  *                                                           iHigh, iHigher)
  117.  *                              double theta, thetacalc();
  118.  */
  119.  
  120. #include <math.h>
  121.  
  122. #if defined(WIN32)
  123. #define PI 3.14159265358979
  124. #endif
  125. #define PID2 1.570796327
  126.  
  127. double
  128. thetacalc (data, iLower, iLow, iHigh, iHigher)
  129.      struct point *data;        /* sequence of coordinates of lines */
  130.      long iLower, iLow,         /* data indices of lower W-segment */
  131.        iHigh, iHigher;          /* data indices of higher W-segment */
  132. {
  133.   long deltaX, deltaY;          /* x,y lengths of segments */
  134.   double gammaLow, gammaHigh,   /* angle of low and high segments */
  135.     angle,                      /* theta angle calculation */
  136.     k,                          /* constant */
  137.     atan ();
  138.  
  139. /* calculate angle of lower segment */
  140.   deltaX = data[iLow].x - data[iLower].x;
  141.   deltaY = data[iLow].y - data[iLower].y;
  142.   if (deltaY >= 0 && deltaX >= 0)
  143.     k = 0.0;
  144.   else if (deltaY >= 0 && deltaX < 0)
  145.     k = PI;
  146.   else if (deltaY < 0 && deltaX >= 0)
  147.     k = 0.0;
  148.   else if (deltaY < 0 && deltaX < 0)
  149.     k = PI;
  150.   if (deltaX == 0)
  151.     gammaLow = (deltaY >= 0) ? PID2 : -(PID2);
  152.   else
  153.     gammaLow = atan (((double) deltaY / (double) deltaX)) + k;
  154.  
  155. /* calculate angle of higher segment */
  156.   deltaX = data[iHigher].x - data[iHigh].x;
  157.   deltaY = data[iHigher].y - data[iHigh].y;
  158.   if (deltaY >= 0 && deltaX >= 0)
  159.     k = 0.0;
  160.   else if (deltaY >= 0 && deltaX < 0)
  161.     k = PI;
  162.   else if (deltaY < 0 && deltaX >= 0)
  163.     k = 0.0;
  164.   else if (deltaY < 0 && deltaX < 0)
  165.     k = PI;
  166.   if (deltaX == 0)
  167.     gammaHigh = (deltaY >= 0) ? PID2 : -(PID2);
  168.   else
  169.     gammaHigh = atan (((double) deltaY / (double) deltaX)) + k;
  170.  
  171.   angle = gammaHigh - gammaLow;
  172.   if (angle > PI)
  173.     angle = angle - 2.0 * PI;
  174.   else if (angle < -PI)
  175.     angle = 2.0 * PI + angle;
  176.  
  177.   return (angle);
  178. }
  179.  
  180.  
  181. /* THETASMTH:   function smooths theta plot with median filter
  182.  *                    usage: thetasmth (theta, nTheta, wwParM)
  183.  *              Note that values of theta at beginning and end are
  184.  *              set to the first and last smoothed values respectively.
  185.  */
  186.  
  187. long
  188. thetasmth (theta, nTheta, wwParM)
  189.      struct theta *theta;       /* theta plot */
  190.      long nTheta;               /* no.values in theta plot */
  191.      long wwParM;               /* ww middle parameter */
  192. {
  193.   double *thetaS;               /* smoothed theta plot */
  194.   long iStrt, iEnd,             /* start and end output theta indices */
  195.     halfWidth,                  /* half width of smooth filter -1 */
  196.     i;
  197.   double median ();
  198.  
  199. /* allocate memory for smoothed theta output */
  200.   if ((thetaS = (double *)
  201.        calloc (nTheta, sizeof (*thetaS))) == NULL) {
  202.     printf ("CALLOC: not enough memory -- sorry");
  203.     return (-2);
  204.   }
  205.  
  206. /* smooth theta plot starting and ending with first correct theta values */
  207.   halfWidth = wwParM / 2;
  208.   for (iStrt = 0; iStrt < nTheta; iStrt++)
  209.     if (theta[iStrt].theta != BORDERTHETA)
  210.       break;
  211.   iStrt = iStrt + halfWidth;
  212.   for (iEnd = nTheta - 1; iEnd >= 0; --iEnd)
  213.     if (theta[iEnd].theta != BORDERTHETA)
  214.       break;
  215.   iEnd = iEnd - halfWidth;
  216.  
  217.   for (i = iStrt; i <= iEnd; i++)
  218.     thetaS[i] = median (&(theta[i - halfWidth]), wwParM);
  219.  
  220.   for (i = iStrt; i <= iEnd; i++)
  221.     theta[i].theta = thetaS[i];
  222.  
  223.   return (0);
  224. }
  225.  
  226.  
  227. /* MEDIAN:      function returns median value of window on theta plot
  228.  *                    usage: m = median (&(theta[i].theta, width)
  229.  *                              double m, median();
  230.  */
  231.  
  232. double
  233. median (theta, width)
  234.      struct theta *theta;       /* windowed portion of theta plot */
  235.      long width;                /* width of smoothing window */
  236. {
  237.   int arrayle[MMAX],            /* array of number of values that is less
  238.                                  * * than or equal to the corresponding
  239.                                  * * <value> */
  240.     widthD2;                    /* window width divided by 2 */
  241.  
  242.   double value[MMAX];           /* value corresponding to less than or
  243.                                  * * equal to array <arrayle> */
  244.  
  245.   register int sumle,           /* running sum of less than or equals */
  246.     sum,                        /* sum of numbers less than or equal to */
  247.     i, j;
  248.  
  249.  
  250. /* initialize */
  251.   for (i = 0; i < MMAX; i++)
  252.     arrayle[i] = 0;
  253.  
  254. /* construct less than or equal to array */
  255.   for (i = 0; i < width; i++) {
  256.     sumle = 0;
  257.     for (j = 0; j < width; j++)
  258.       if (theta[i].theta <= theta[j].theta)
  259.         sumle++;
  260.     arrayle[sumle]++;
  261.     value[sumle] = theta[i].theta;
  262.   }
  263.  
  264. /* find median */
  265.   widthD2 = width / 2;
  266.   sum = 0;
  267.   for (i = 1; i <= width; i++) {
  268.     sum += arrayle[i];
  269.     if (sum > widthD2)
  270.       break;
  271.   }
  272.  
  273.   return (value[i]);
  274. }
  275.  
  276.  
  277. /* WWENDS:      function finds points on curve of equal arc length
  278.  *            above and below a given curve point
  279.  *                  usage: returnFlag = wwends (data, edge, theta, iMidPt,
  280.  *                                                  length, &iLow, &iHigh)
  281.  *              <returnFlag> = -1 if lower end of data reached; = -2
  282.  *              if higher end of data reached; = -3 if both reached;
  283.  *              = 0 if neither reached
  284.  */
  285.  
  286. #include <images.h>             /* image format information */
  287.  
  288. long
  289. wwends (data, edge, theta, iMidPt, length, iLow, iHigh)
  290.      struct point *data;        /* sequence of coordinates of curve */
  291.      struct edge edge;          /* strt and end indices of edges in data */
  292.      struct theta *theta;       /* theta plot curv. and arc length array */
  293.      long iMidPt,               /* middle point of arc */
  294.        length,                  /* arc length from lower to higher point */
  295.       *iLow, *iHigh;            /* lower, and higher indices of data pts */
  296. {
  297.   double halfLength,            /* half length of line around point */
  298.     arcLength,                  /* cumulative length of line */
  299.     arcLengthB4;                /* cumul. length before latest increment */
  300.   long returnFlag,              /* indicates either end of data reached */
  301.     j;
  302.  
  303.   halfLength = (double) length / 2.0;
  304.  
  305. /* find lower point */
  306.   j = iMidPt;
  307.   returnFlag = 0;
  308.   arcLength = 0.0;
  309.   do {
  310.     --j;
  311.     if (j < edge.iLow) {
  312.       j = edge.iLow;
  313.       returnFlag = -1;
  314.       break;
  315.     }
  316.     arcLengthB4 = arcLength;
  317.     arcLength = theta[iMidPt - edge.iLow].length
  318.       - theta[j - edge.iLow].length;
  319.   } while (arcLength < halfLength);
  320.   *iLow = ((arcLength - halfLength) > (halfLength - arcLengthB4))
  321.     ? (++j) : j;
  322.  
  323. /* find higher point */
  324.   j = iMidPt;
  325.   arcLength = 0.0;
  326.   do {
  327.     j++;
  328.     if (j > edge.iHigh) {
  329.       j = edge.iHigh;
  330.       returnFlag += -2;
  331.       break;
  332.     }
  333.     arcLengthB4 = arcLength;
  334.     arcLength = theta[j - edge.iLow].length
  335.       - theta[iMidPt - edge.iLow].length;
  336.   } while (arcLength < halfLength);
  337.   *iHigh = ((arcLength - halfLength) > (halfLength - arcLengthB4))
  338.     ? (--j) : j;
  339.  
  340.   return (returnFlag);
  341. }
  342.